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Abstract 

A fast algorithm for B-splines in mixed models is presented. B-splines have 
local support and are computational attractive, because the corresponding 
matrices are sparse. A key element of the new algorithm is that the local 
character of B-splines is preserved, while in other existing methods this local 
character is lost. The computation time for the fast algorithm is linear in 
the number of B-splines, while computation time scales cubically for existing 
transformations. 
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1 Introduction 


Penalized regression using B-splines can be computationally efficient, be¬ 
cause of the local character of B-splines. The corresponding linear equations 
are sparse and can be solved quickly. However, the main problem is to fold 
the optimal value for the penalty parameter. A good way to approach this 
problem is to use mixed models and restricted maximum likelihood (REML; 
Patterson & Thompson, 1971). Several methods have been proposed to 


transform the original penalized B-spline model to a mixed model (Currie & 


Durban 2002 Lee & Durban, 2011). A problem with existing transforma¬ 


tions to mixed models is that the local character of the B-splines is lost, which 
reduces the computational efficiency. For relatively small datasets this is not 
a major issue. However, for long time series, for example with measurements 
every five minutes for several months, the computational efficiency becomes 
quite important. 

In this paper I present a new transformation to a mixed model. This model 
is closely related to the transformation proposed by Currie & Durban (2002). 


However, the computation time in the transformation of Currie & Durban 


(2002) increases cubically in the number of B-splines, while for the new trans¬ 
formation the computation time increases linearly in time, using sparse ma¬ 
trix algebra (Furrer & Sain, 2010). One of the key elements of the pro¬ 


posed algorithm is that the transformation preserves the local character of 
B-splines, and all the equations can be solved quickly. 

The paper is organized as follows. In Section [2] relevant information about B- 
splines is given. In Section ([3]) first the P-splinc model (Eilers & Marx, 1996) 
is described and the transformation to mixed models by Currie & Durban 


(2002) is stated. The new transformation is presented, and details for a 
sparse mixed model formulation are given. In Section [4] the R-code is briefly 
described and a comparison is made between the computation time of the 
new method and the transformation by Currie & Durban (2002). 


2 B-splines 


In this preliminary section, a few relevant details about B-splines are given. 


For a detailed overview of B-splines, see for example De Boor (1978) and 

B-splines have local support 


Hastie et al. (2009) 


tant and can speedup calculations considerably. 


This property is irnpor- 
To illustrate the idea of 
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Figure 1: Quadratic B-spline basis for the interval [x m i n ,x max ] = [0,10], formed 
by B 1,2(3;), B 2,2(2;), ... , -Bi2,2(.t). The distance between the splines is unity. The 
red curves are second-order derivatives of quartic B-splines; for clarity only the 
first and last one are shown. The second-order derivatives of quartic B-splines can 
be constructed from quadratic B-splines, see equation (fll). 


the local support, see Figure 1, with quadratic (i.e. degree q = 2) B-splines. 
Throughout the paper we will assume equal distance between the splines, de¬ 
noted by h. I 11 the example presented in Figure 1 the distance is unity, h — 1. 
The number of B-splines will be denoted by m. The first quadratic B-spline, 
B i )2 (t) is zero outside the interval [—2,1], The last one, B 12,2 (t), is zero 
outside the interval [9,12], So, for this example, there are m = 12 quadratic 
B-splines which define the B-spline basis for the domain [x m i n , x mSLX ] = [0,10] 
of interest. 


The second derivative of B-splines is given by (De Boor, 1978): 


h2B "j, q ( x ) = B j,g- 2 O) - 2 B j+hq _ 2 (x) + B j+ 2 ,,_ 2 (x) . (1) 


The second derivative is illustrated in Figure 1. The red curves are second- 
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order derivatives of quartic (fourth-degree) B-splines. As can be seen from 
equation (|TJ) it can represented as a linear combination of three quadratic 
B-splines. 


3 P-splines and Mixed Models 

In this section we will give a brief description of P-splines (Eilcrs & MarxJ 


1996). Let n be the number of observations. Suppose the variable y = 


(j/i,..., y n )' depends smoothly on the variable x = (aq,..., x n )'. Let B = 
(Bi > 2 (x), ..., -B mj2 (x)) be a nx m matrix, and a = (aq, a 2 ,..., a m )' be a 
vector of regression coefficients. Then the following objective function to be 
minimized can be defined: 

5(a) = (y - Ba)'(y - Ba) + A a'D'Da , (2) 

where A > 0 is a penalty or regularization parameter, and D is an (m — 2) x m 
second-order difference matrix (see e.g. Eilers & Marx, 1996). 

Currie & Durban (2002) showed that equation (|2j) can be reformulated as a 


mixed model: 


y = Xb + Zu + e , u ^ N( 0, —Q V 2 ) , e ^ N( 0.1 a 

A 


(3) 


where X and Z are design matrices, Q is a precision matrix, b = (6q, bi)' are 
the fixed effects, u = (iq, w 2 ,..., w m _ 2 )' the random effects, e is the residual 
error, and a 2 is the residual variance. 


Currie & Durban (2002) used the following transformation: 

a = Gb + D'(DD'' _1 ’ 


u 


(4) 


where G is an m x 2 matrix with columns go = (1,1,..., 1)' and gi = 
(1, 2,..., m)'. This transformation gives the following expressions for the 
design matrices and the precision matrix: 


X = BG , Z = BD'(DD 


A-l 


Q = I 


(5) 


The mixed model equations (Henderson, 1963) corresponding to equation ([3]) 
are given by: 

'x'x xz A fb\ = /xy 

Z'X Z'Z + AqMuJ Iz'y 
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( 6 ) 





















( 7 ) 


The coefficient matrix Ca in equation (J6]) is given by: 


/X'X XZ \ 

^Z'X Z'Z + AQy ' 


This coefficient matrix is dense, since the local character of the B-splines has 
been destroyed by equation (J5]) . This implies that the computation complex¬ 
ity for solving equation (J6| is (9(m 3 ). 

The following transformation preserves the local character of the B-splines: 


a = Gb + D 7 u . 


( 8 ) 


Figure 1 illustrates the underlying idea of this transformation. The quadratic 
B-splines basis consists of m = 12 B-splines. This quadratic B-spline basis is 
transformed to a second-order derivative quartic B-splines basis of m — 2 = 
10 B-splines, plus a parameter for intercept and linear trend b\. The 
second-order derivative quartic B-splines can be constructed from quadratic 
B-splines by second-order differencing. Using the new transformation the 
design and precision matrices are given by: 


X = BG , Z = BD' , Q = DD'DD' . 


(9) 


Let us refer to equations ([3]) and 0 as a Mixed Model of B-splines (MMB), 
since it uses the B-splines directly as building blocks for the mixed model. 
The matrix Z'Z + AQ has bandwidth 4. This implies that Ca is sparse and 
computation complexity has been reduced to 0{m). An efficient way to 


calculate the REML profile log likelihood (Gilmour et al., 1995 Crainiceanu 


& Ruppert, 2004 Searle et al. 2009) is given by the following four steps : 


1. Sparse Cholesky factorization (Furrer & Sain, 2010): 


C A = U a U a , 


where Ua is an upper-triagonal matrix. 


2. Forward-solve and back-solve (Furrer & Sain 2010), with w a vector 
of length m : 


U A w = 


X'y 

Z'y 


U, 


u 


= w . 


( 10 ) 


3. Calculate a 2 (Johnson & Thompson, 1995), p is the dimension of the 
fixed effects: 


a 2 = 


y'y - b'X'y - u'Z'y 

n — p 


( 11 ) 
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4. 

REML log profile likelihood ( 

Gilmour et al. 

1995; 

Crainiceanu & Rup- 


pert, 

2004; 

Searle et ah, 

2009): 


L(X) = —^ (2 log |Ua| — (m — p) log X + (n — p) loga 2 + C) , (12) 


where C is a constant: C = n — p — log |Q|. 

A one-dimensional optimization algorithm can be used to find the maximum 
for L(X). The computation time is linear in m. 


4 R-package MMBsplines 


An R-package, MMBsplines, is available at GitHub: 
https://github.cora/martinboer/MMBsplines.git, 

The sparse matrix calculations are done with the spam package Furrer & 
Sain (2010). The B-splines are constructed with splineDesignO of the 


splines library. 

The following example code sets some parameter values and runs the simu¬ 
lations: 


nobs = 1000; xmin = 0; xmax = 10 
set.seed(949030) 

sim.fun = function(x) { return(3.0 + 0.1*x + sin(2*pi*x))} 
x = runif(nobs, min = xmin, max = xmax) 
y = sim.fun(x) + 0.5*rnorm(nobs) 

A fit to the data on a small grid can be obtained as follows, using m = 100 
quadratic B-splines: 

obj = MMBsplines(x, y, xmin, xmax, nseg = 100) 

xO = seqCxmin, xmax, by=0.01) 

yhat = predict(obj, xO) 

ylin = predict(obj, xO, linear = TRUE) 

ysim = sim.fun(xO) 

Figure [2] shows the result, with A max = 1.33. 

The Currie and Durban transformation can be run by setting the sparse argu¬ 
ment to FALSE: 

obj = MMBsplines(x, y, xmin, xmax, nseg = 100, sparse = FALSE) 
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Figure 2: Fit of the simulated data using MMBsplines, with A max = 1.33. The 
blue line is the true simulated line, the red line is the fitted value. The green line 
is the linear trend. 


For m = 100, as in Figure [2j the differences in computation time are small. If we 
increase the length of the simulated time series, with a fixed stepsize h = 0.1, the 
advantage of the MMB-splines method becomes clear, see Figure 3. As expected 
the Currie and Durban transformation computation time increases cubical in the 
number of B-splines, computatation time for MMB-splines is linear in m. 


5 Conclusion 


The MMB-splines method presented in this paper seems to be an attractive way 
to use B-splines in mixed models. The method was only presented for quadratic 
splines, but also cubical or higher-degree B-splines could have been used. Other 


generalization s are a lso possible, for example extension to multiple penalties (Cur- 
rie & Durban, 2002) or multiple dimensions (Rodriguez-Alvarez et al., 2014). 
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Figure 3: Comparison of computation times. The computation time for the Currie 
and Durban transformation is cubical in the number of B-splines. The computation 
time for MMB-splines is linear in the number of B-splines. 
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